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Abstract 

In many signal processing applications, one wishes to acquire images that are sparse in trans- 
form domains such as spatial finite differences or wavelets using frequency domain samples. For such 
£SJ applications, overwhelming empirical evidence suggests that superior image reconstruction can be 

t-H obtained through variable density sampling strategies that concentrate on lower frequencies. The 

wavelet and Fourier transform domains are not incoherent because low-order wavelets and low-order 
frequencies are correlated, so compressed sensing theory does not immediately imply sampling strate- 
gies and reconstruction guarantees. In this paper we turn to a more refined notion of coherence - 
the so-called local coherence - measuring for each sensing vector separately how correlated it is to 
the sparsity basis. For Fourier measurements and Haar wavelet sparsity, the local coherence can be 
controlled, so for matrices comprised of frequencies sampled from suitable power-law densities, we can 
prove the restricted isometry property with near-optimal embedding dimensions. Consequently, the 
variable-density sampling strategies we provide — which are independent of the ambient dimension 
i | up to logarithmic factors — allow for image reconstructions that are stable to sparsity defects and 

robust to measurement noise. Our results cover both reconstruction by ^-minimization and by total 
variation minimization. 

, ^, 1 Introduction 

The measurement process in a wide range of imaging applications such as radar, sonar, astronomy, and 
computer tomography, can be modeled - after appropriate approximation and discretization - as taking 
samples from weighted discrete Fourier transforms |T6] . Similarly, it is well known in the medical imaging 
literature that the measurements taken in Magnetic Resonance Imaging (MM) are well modeled as Fourier 
coefficients of the desired image. Within all of these scenarios, one seeks strategies for taking frequency 
domain measurements so as to reduce the number of measurements without degrading the quality of 
I image reconstruction. A central feature of natural images that can be exploited in this process is that 

£SJ they allow for approximately sparse representation in suitable bases or dictionaries. 

The theory of compressed sensing, also termed compressive sampling, as introduced in |14[ [8] , fits 
comfortably into this set-up: its key observation is that signals which allow for a sparse or approximately 
sparse representation can be recovered from relatively few linear measurements via convex approximation, 
provided these measurements are sufficiently incoherent with the basis in which the signal is sparse. 
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1.1 Imaging with partial frequency measurements 

Much work in compressed sensing has focused on the setting of imaging with frequency domain mea- 
surements |Hl 13 [34] and, in particular, towards accelerating the MRI measurement process (221 [21]. For 
images that are localized in the pixel domain, the incoherence between pixel and Fourier bases implies 
that uniformly subsampled discrete Fourier transform measurements can be used to achieve near-optimal 
oracle reconstruction bounds: up to logarithmic factors in the discretization size, any image can be ap- 
proximated from s such frequency measurements up to the error that would be incurred if the image 
were first acquired in full, and then compressed by setting all but the s largest-magnitude pixels to zero 
[3U [311 [71 [9]. Compressed sensing recovery guarantees hold more generally subject to incoherence be- 
tween sampling and sparsity transform domains. Unfortunately, natural images are generally not directly 
sparse in the pixel basis, but rather with respect to transform domains more closely resembling wavelet 
bases. As low-scale wavelets are highly correlated (coherent) with low frequencies, sampling theorems for 
compressive imaging with partial Fourier transform measurements have remained elusive. 
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A number of empirical studies, including the very first papers on compressed sensing MRI [2"2l 121) . 
suggest that better image restoration is possible by subsampling frequency measurements from variable 
densities preferring low frequencies to high frequencies. In fact, variable-density MRI sampling had been 
proposed previously in a number of works outside the context of compressed sensing, although there 
did not seem to be a consensus on an optimal density [23] |36l [38] [2Sj HB] ■ For a thorough empirical 
comparison of different densities from the compressed sensing perspective, see [55] , 

Guided by these observations, the authors of [50] estimated the coherence between each element of 
the sensing basis with elements of the sparsity basis separately as a means to derive optimal sampling 
strategies in the context of compressed sensing MRI. In particular, they observe that incoherence-based 
results in compressed sensing imply exact recovery results for more general systems if one samples a row 
from the measurement basis proportionally to its squared maximal correlation with the sparsity-inducing 
basis. For given problem dimensions, they find the optimal distribution as the solution to a convex 
problem. In |29) . the approach of optimizing the sampling distribution is combined with premodulation 
by a chirp, which by itself is another measure to reduce the coherence [5] [30]. A similar variable-density 
analysis already appeared in |33j in the context of sampling strategies and reconstruction guarantees 
for functions with sparse orthogonal polynomial expansions and will also be the guiding strategy of this 
paper (cf. Section [jj]). 

1.2 Contributions of this paper 

In this paper we derive near-optimal reconstruction bounds for variable-density subsampled discrete 
Fourier transform measurements for both wavelet sparsity and gradient sparsity models. More precisely, 
up to logarithmic factors in the discretization size, any image can be approximated from s such mea- 
surements up to the error that would be incurred if the wavelet transform or the gradient of the image, 
respectively, were first computed in full, and then compressed by setting all but the s largest-magnitude 
coefficients to zero. Note that the reconstruction results that have been derived for uniformly-subsamplcd 
frequency measurements [5] [TS] are not optimal in this sense. 

A major role in determining the best sampling density will be played by the local coherence of the 
sensing basis with respect to the sparsity basis, as introduced in Section [5j Consequently, an important 
ingredient of our analysis is Theorem |6.2[ which provides frequency-dependent bounds on inner products 
between rows of the orthonormal discrete Fourier transform and rows of the orthonormal discrete Haar 
wavelet transform. In particular, the maximal correlation between a fixed row in the discrete Fourier 
transform and any row of the discrete Haar wavelet transform decreases according to an inverse power 
law of the frequency, and decays sufficiently quickly that the sum of squared maximal correlations scales 
only logarithmically with the discretization size N. This implies, according to the techniques used in 
[551 1521 0] , that subsampling rows of the discrete Fourier matrix proportionally to the squared correlation 
results in a matrix that has the restricted isometry property of near-optimal order subject to appropriate 
rescaling of the rows. 

For reconstruction, total variation minimization [35] [5] [27] [32 [Til HO] will be our algorithm of choice. 
In the papers |25j and [26], total variation minimization was shown to provide stable and robust image 
reconstruction provided that the sensing matrix is incoherent with the Haar wavelet basis. Following the 
approach of |25j . we prove that from variable density frequency samples, total variation minimization can 
be used for stable image recovery guarantees. Due to the effects of variable density sampling we will no 
longer need the samples to be incoherent to the Haar basis. 

1.3 Outline 

The remainder of this paper is organized as follows. Preliminary notation is introduced in Section 2. 
The main results of this paper along with a numerical illustration are contained in Section 3. Section 
4 reviews compressed sensing theory and Section 5 presents recent results on sampling strategies for 
coherent systems. The main results on the coherence between Fourier and Haar wavelet bases is provided 
in Section 6, and proofs of the main results are contained in Section 7. We conclude with a summary 
and a discussion of open problems in Section 8. 
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2 Preliminaries 



2.1 Notation 

In this paper, we consider discrete images, that is, TV x A blocks of pixels, and represent them as discrete 
functions / £ <C NxN . We write /(fi,f 2 ) to denote any particular pixel value and for a positive integer 
A, we denote the set {1,2,..., A} by [A]. By fi o f 2 we denote the Hadamard product, i.e., the image 
resulting from pointwise products of the pixel values, /i o /2(ii,<2) = ^2)^2 (^i , £2)- On the space 

of such images, the £ p vector norm is given by ||/|| p = (J2t ± t 2 l/(*i)*2)| p ) 1 ^ P , 1 < p < 00, and ||/||oo = 
max( tli(2 ) |/(fi,f 2 )|- The inner product inducing the £ 2 vector norm is (f,g) — Y]f r t2 f i^i ^'2)9(^1^2) , 
where z denotes the complex conjugate of number z £ C. By an abuse of notation, the "i?o-norm" 
ll/llo = #{(ii)*2) : /(*i) ^2) 0} counts the number of non-zero entries of /. 

An image / is called s-sparse if ||/|| < s. The error of best s-term approximation of an image / in 
l p is defined as 

°"*(/)p = n in f 11/ "Slip- 

Clearly, a s (f) p = if / is s-sparse. Informally, / is called compressible if er s (/)i decays quickly as s 
increases. 

For two nonncgative functions /(f) and g{t) on the real line, we write / > g (or / < g ) if there exists 
a constant C > such that f(t) > Cg(t) (or /(f) < Cg(t), respectively) for all f > 0. 
The discrete directional derivatives of / £ C NxN are defined pixel- wise as 



f x eC N - lxN , / x (fi,f 2 ) = /(fi + l,f 2 ) -/(fi, f 2 ) 
f y eC NxN -\ fyihM) = f(ti,t 2 + l)-f(t 1 ,t 2 ) 
The discrete gradient transform V : C NxN — > C WxArx2 i s defined in terms of the directional derivatives 



via 



V/(fi,f 2 ) := (/x(tl,t 2 ), fy{tl,t 2 )), 



where the directional derivatives are extended to x by adding zero entries. The fofa/ variation 
semi-norm is the t\ norm of the image gradient, 

\\K\tv :=||V/|| 1 = ^(|/,(f 1 ,f 2 )| + |/ ?/ (f 1 ,f2)|). 

tl,*2 

Here we note that our definition is the anisotropic version of the total variation semi-norm. The 
isotropic total variation semi-norm becomes the sum of terms 

|/x(*i,*2) + */„(*i,* 2 )| - (fxihMf + fyihM) 2 ) 1 ' 2 - 

The isotropic and anisotropic total variation semi-norms are thus equivalent up to a factor of y/2. 



2.2 Bases for sparse representation and measurements 



The Haar wavelet basis is a simple basis which allows for good sparse approximations of natural images. 
We will work primarily in two dimensions, but first introduce the univariate Haar wavelet basis as it will 
nevertheless serve as a building block for higher dimensional bases. 

Definition 2.1 (Univariate Haar wavelet basis). The univariate discrete Haar wavelet system is an 
orthonormal basis of C 2 " consisting of the constant function h°(t) = 2~ p / 2 , the step function /ij = h 1 
given by 

h i M _j 2-P/ 2 , l<f<2f-\ 
{ } ~ \ -2-P/ 2 , 2P- 1 <t<2P, 

along with the dyadic step functions 



2t/i 1 (2"f - 


-l) 






for 


i2 P-n < t <• £ 2 P- n + 2P-™- 1 


< -2^ 


for 


£2P-n + 2 P-n-l < t < £ 2 P- n + 2 p ~ n 




V 


else, 





for (n, I) £ I? satisfying < n < p and < I < 2". 
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To define the bivariate Haar wavelet basis of C 2Px2P , we extend the univariate system by the window 
functions 



h° n At) = 2%h°(2 n t-l) = 



2 V for £2P- n <t<(£+ l)2P- n 

else. 



The bivariate Haar wavelet system can now be defined via tensor products of functions in the extended 
univariate system. In order for the system to form an orthonormal basis of C 2 x2 , only tensor products 
of univariate functions with the same scaling parameter n are included. 

Definition 2.2 (Bivariate Haar wavelet basis). The bivariate Haar system of C 2Px2P consists of the 
constant function /i^ 0,0 ) given by 

h^(t ll t 2 ) = h (t 1 )h°(t 2 ) = 2-P 

and the functions h e n e with indices in the range 0<n<p,f=(fi,<2)€Z 2 n [0, 2™) 2 , and 
e = (ei,e 2 ) E {{0, 1}', {1,0}, {1, 1}} given by 

K 4 {UM) = K\ li {t l )hi t2 {t 2 ). 

We denote by % the bivariate Haar transform / — > ( (/, h n j) ) n( and, by a slight abuse of notation, also 
the unitary matrix representing this linear map. 

We will also work with discrete Fourier measurements. 

Definition 2.3 (Discrete Fourier basis). Let N — 2 P . The one-dimensional discrete Fourier system is an 
orthonormal basis of C N consisting of the vectors 

<p k (t) = -^=e i2ntk/N , -N/2 + l<t<N, 



indexed by discrete frequencies in the range —N/2 + 1 < k < N/2. The two-dimensional discrete Fourier 
basis of C NxN [ s j us t a tensor product of one-dimensional bases, namely 

Vk lM (t 1 M) = ^e a < t ^+ t ^l N 1 -N/2 + 1 <t l7 t 2 < N/2, (2.3) 

indexed by discrete frequencies in the range —N/2 + 1 < k\, k 2 < N/2. 

We denote by T the two-dimensional discrete Fourier transform / — > ( (/, tp^ ,fc 2 ) ) fel k2 an d, again, also 
the associated unitary matrix. Finally, we denote by Tq its restriction to a set of frequencies ft C [N] 2 . 



3 Main results 

Our main results say that appropriate variable density subsampling of the discrete Fourier transform will 
with high probability result in a set of measurements admitting stable image reconstruction via total 
variation minimization or ^-minimization. 

While our recovery guarantees are robust to measurement noise, the model we consider corresponds 
to a weighted ^2-norm such that high-frequency measurements have higher sensitivity to noise. Loosely 
speaking, this compensates for the smaller number of samples in regions of low sampling probability. Our 
first result concerns stable recovery guarantees for total variation minimization. 

Theorem 3.1. Fix integers N = 2 p 1 m 1 and s such that s > log(iV) and 

m > s log 3 (s) log 5 (TV). 
Select m frequencies {(lu : ' 1i ui 2 ')}"L 1 C {— N/2 + 1, . . . , 7V/2} 2 i.i.d. according to 

Prob[(wj,wg) = (A;i,fc2)] = C N min (c, ^^2 ) =-V{h,k 2 ), -N/2 + 1 < fci, k 2 < N/2, (3.1) 

where C is an absolute constant and CV is chosen such that r\ is a 'probability distribution. 

Consider the weight vector p = (pj)JL 1 with pj = (l/rj(uj{, o^)) 1 / 2 , and assume that the noise vector 
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£ = satisfies \\p o £|| 2 < E\fm, for some e > 0. TTien luif/i probability exceeding I — N clo s ( s ) j 

t/ie following holds for all images f € C NxN : 

Given noisy partial Fourier measurements y = J 7 ^/ + £ ; ^ e estimation 

f* = argmin ||ff||rv such that \\p o (J^ n g - y)\\ 2 < e-s/m, (3.2) 

approximates f up to the noise level and best s-term approximation error of its gradient: 

V-f#b< g^gMi + . 



Disregarding measurement noise, the error rate provided in Theorem 3.1 (and also the one in Theo- 
rem 3.2 below) is optimal up to logarithmic factors in the ambient image dimension. This follows from 
classical results about the Gel'fand width of the £i-ball due to Kashin [TS] and Garnaev-Gluskin [T7] , 

Numerical results such as those detailed in [35] and illustrated below in Figure [I] confirm that variable- 
density sampling strategies significantly outperform uniform sampling strategies as well as deterministic 
sampling strategies, and Theorem |3 . 1 1 provides theoretical justification for such observations. 

Our second result focuses on stable image reconstruction by ^-minimization in the Haar wavelet 
transform domain. It is a direct consequence of applying the Fourier- wavelet incoherence estimates 



derived in Theorem 15.21 to Theorem 16.21 

Theorem 3.2. Fix integers N = 2 p ,m, and s such that s > log (AT) and 

m > s log 2 (AO log 3 (s). 

Select m frequencies f2 = {(cj^, w^jTLi <- {~ A^/2 + 1, . . . , Af/2} 2 i.i.d. according to the density rj as in 



(3.1 1 and assume again that the noise vector £ = (£j)jLi satisfies the weighted l^-constraint with weight p 

Then with probability exceeding 1 — _/V~ c ' los3 ( s ) ; the following holds 



3.1 



and noise level e as in Theorem 

for all images f £ C NxN : Given noisy measurements y = J-q/ + £, the estimation 



f* = argmin ||W<?||i such that \\po (Tng — y)\\2 < £ 
approximates f up to the noise level and best s-term approximation error in the bivariate Haar basis: 

W- (uf) s \\i , _ 



l/-/ # ll2< 



Even though the required number of samples m in Theorem |3.2| is smaller than the number of samples 
required for the total variation minimization guarantees in Theorem 3.1 we find that total variation 
minimization requires fewer measurements empirically. This may be due to the fact that the gradient of 
a natural image has stronger sparsity than its Haar wavelet representation. For this reason we focus on 
total variation minimization. Independent of this observation, we strongly suspect that the additional 
logarithmic factors in the number of measurements stated in Theorem |3.1| are an artifact of the proof, 
and that it should be possible to strengthen the result to obtain a similar recovery guarantee with the 
number of measurements as in Theorem [3^21 



4 Compressed sensing background 
4.1 The restricted isometry property 

Under very mild assumptions on the matrix $ : C N — > C m , any fc-sparse x € can be recovered from 
y = $a; as the solution to the optimization problem: 

x = argmin |jz||o such that $z = y 

One of the fundamental results in compressed sensing is that this optimization problem, which is NP-hard 
in general, can be relaxed to an ^-minimization problem if one asks that the matrix $ restricted to any 
subset of 2k columns be well-conditioned. This property is quantified via the so-called restricted isometry 
property as introduced in [3] : 
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Definition 4.1 (Restricted isometry property). Let $ £ C nxN . For s < N, the restricted isometry 
constant 8 S associated to <& is the smallest number 5 for which 

(l-6)\\x\\l<\\$x\\i<(l + d)\\x\\ 2 2 

for all s-sparse vectors x £ C N . If 5 S < 5, one says that $ has the restricted isometry property (RIP) of 
order s and level 5. 

The restricted isometry property ensures stability: not only sparse vectors, but also compressible 
vectors can be recovered from the measurements via ^-minimization. It also ensures robustness to 
measurement errors. 

Proposition 4.2 (Sparse recovery for RIP matrices). Assume that the restricted isometry constant 5§ s 
of $ £ C mx,v satisfies 85s < |. Let x £ C N and assume noisy measurements y = &x + £ with ||£||2 < £■ 
Then 

x^ — arg min ||^||i subject to \\<&z — j/| 1 2 < £ 

satisfies 

k - a;* 2 < W— + £• 

In particular, reconstruction is exact, x& — x, if x is s-sparse and e = 0. 

There are stronger versions of this result which allow for weaker constraints on the restricted isometry 
constant |24j . However, our version is a corollary of the following proposition, which generalizes the 
results from |8j and seem to have appeared first in [25] . This proposition will also play an important role 
in the proof of our main results. 

Proposition 4.3 (Stable recovery for RIP matrices, [IS]). Suppose that 7 > 1 and $ € i^ mxN satisfies 
the restricted isometry property of order bk'y 2 and level 5 < 1/3, and suppose that u € C N satisfies a tube 
constraint 

\\®uh ~ £ - 

Suppose further that for a subset S of cardinality \S\ — k, the signal u satisfies a cone constraint 

\WsA\i < 7ll«sf||i +£• 



Then 



HI2 



< 



7^ 



(4.1) 



Proposition 



4.2 follows from Proposition 4.3 by noting that the minimality of x# implies a cone con- 
straint for the residual x — x& over the support of the s largest-magnitude entries of x. For completeness, 
the proof of Proposition |4.3| is recalled in the appendix. 



4.2 Bounded orthonormal systems 

While the strongest known results on the restricted isometry property concern random matrices with 
independent entries such as Gaussian or Bernoulli, a scenario that has proven particularly useful for 
applications is that of structured random matrices with rows chosen from a basis incoherent to the basis 
inducing sparsity (see below for a detailed discussion on the concept of incoherence). The resulting 
sampling schemes correspond to bounded orthonormal systems, and such systems have been extensively 
studied in the compressed sensing literature (see |31j for an expository article including many references) . 

Definition 4.4 (Bounded orthonormal system). Consider a set T equipped with probability measure v. 

• A set of functions {ipj : T — > C, j G [N]} is called an orthonormal system with respect to v if 
J T i[)j{x)'4>k(x)dv(x) — Sjk, where Sjk denotes the Kronecker delta. 

• An orthonormal system is said to be bounded with bound K if sup^gnyi ||V'j( a ')||oo < K. 



G 



For example, the basis of complex exponentials ipj( x ) — Gxp (i2njx) forms a bounded orthonormal 
system with optimally small constant K = 1 with respect to the uniform measure on T = {0, jj, . . . , ^r"}, 
and d-dimensional tensor products of complex exponentials form bounded orthonormal systems with 
respect to the uniform measure on the set T d . A random sample of an orthonormal system is the vector 
(ipi(x), . . . , ip]sr(x)), where a; is a random variable drawn according to the associated distribution v. Any 
matrix whose rows are independent random samples of a bounded orthonormal system, such as the 
uniformly subsampled discrete Fourier matrix, will have the restricted isometry property: 

Proposition 4.5 (RIP for bounded orthonormal systems, [3T]). Consider the matrix \t € C mxN whose 
rows are independent random samples of a orthonormal system {ipj, j € [N]} with bound K > 1 and 
orthogonalization measure v. If 

to > <r 2 iY 2 .slog 3 (s)log(iV), 

for some s > log(ivj^J then with probability at least 1 — J\f- C1 °s ( s ) j the restricted isometry constant S s 
of satisfies 5 S < 5. 

An important special case of a bounded orthonormal system arises by sampling a function with sparse 
representation in one basis using measurements from a different, incoherent basis. The mutual coherence 
between a unitary matrix A € i^nxn rows {a,j)^ =l and a unitary matrix B G tr^NxN w ^h rows 

(bj)f =1 is given by 

fi(A,B) = sup | (dj,bk) | 
j,k 

The smallest possible mutual coherence is /i = TV -1 / 2 , as realized by the discrete Fourier matrix and 
the identity matrix. We call two orthonormal bases A and B mutually incoherent if fj, = 0(N^ 1 ^ 2 ) or 
fj, = 0(\og a (N)N~ 1/2 ). In this case, the rows (bj)f =1 of the basis B = \/~NBA* constitute a bounded 
orthonormal system with respect to the uniform measure. Propositions |4.5| and |4.2| then imply that, 
with high probability, signals / £ C N of the form / = Ax for x sparse can be stably reconstructed from 
uniformly subsampled measurements y = Bf = Bx, as B has the restricted isometry property. 

Corollary 4.6 (RIP for incoherent systems, [H]). Consider orthonormal bases A,B£ <£ NxN w ith mutual 
coherence bounded by /i(A,B) < KN~ X / 2 . Fix 8 > and integers N,m, and s such that s > log(A) and 

m > 8- 2 K 2 s log 3 (s)\og(N). (4.2) 

Consider the matrix <J> £ (£ mxN formed by uniformly subsampling m rows i.i.d. from the the matrix 
B = \/NBA*. Then with probability at least 1 — N~ ci ° e W , the restricted isometry constant S s of ^7= < i ) 
satisfies S s < 5. 



5 Local coherence 



The coherence-based sparse recovery results implied by Corollary 4.6 do not take advantage of the 
full range of applicability of bounded orthonormal systems. As argued in (33], Proposition 4.5 implies 



comparable sparse recovery guarantees for a much wider class of sampling/sparsity bases through pre- 
conditioning resampled systems. In the following, we formalize this approach through the notion of local 
coherence. 

Definition 5.1 (Local coherence). The local coherence of an orthonormal basis {fj}^ =1 of C N with 
respect to the orthonormal basis {^fclfcL 
by 



of C is the function /i ($, ty) £ M. defined coordinate-wise 



sup \(tpj,i>k)\- 

Kk<N 



The following result shows that we can reduce the number of measurements to in (4.6 ) by replacing the 
bound K on the coherence in (4.6 ) by a bound on the ^-norm of the local coherence, provided we sample 



rows from $ appropriately using the local coherence function. It can be seen as a direct finite-dimensional 
analog to Theorem 2.1 in [33], but for completeness, we include a short self-contained proof. 



For matrices consisting of uniformly subsampled rows of the discrete Fourier matrix, it has been shown in |12| that this 
constraint is not necessary. 
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Theorem 5.2. Let <£> = {<^j}jLi an d ^ — {"0fc}feLi be orthonormal bases of C N . Assume the local 
coherence of <£> with respect to "J is pointwise bounded by the function k, that is sup \(ipj,ipk}\ < Kj. 

l<k<N 

Let s > log(iV) ; suppose 

m> S- 2 \\k\\ 2 2 s log 3 (s)\og(N), 

and choose m (possibly not distinct) indices j G O C [iV] z.icL /rom i/ie probability measure v on [N] 
given by 

Consider the matrix A € C mxN with entries 

A jtk = (<Pj,ipk), jen,ke[N], 

and consider the diagonal matrix D = diag(d) G C N with dj = Wk^/kj. Then with probability at least 
1 — N~ clog3l - 8 \ the restricted isometry constant S s of the preconditioned matrix -j=DA satisfies S s < S. 

Proof. We show that the system {(fj} — {djtpj} is an orthonormal system with respect to v in the sense 
of Definition |4.4| Indeed, 

N ff ii ii ii ii 2 

' \K\\2 A /||k||2 ^\ K j 



3=1 J J 

iY 

£^(fci)^(fc 2 ) = <S fcl ,fc 2 



3=1 j=i 



3=1 

hence the ^ form an orthonormal system with respect to v. Noting that this system is bounded with 



bound \\k U, the result follows from Proposition 4.5 □ 



Remark 5.3. Note that the local coherence not only appears in the embedding dimension m, but also 
in the sampling measure. Hence a priori, one cannot guarantee the optimal embedding dimension if one 
only has suboptimal bounds for the local coherence. That is why the sampling measure in Theorem |5.2| 
is defined via the (known) upper bounds k and ||/c||2 rather than the (usually unknown) exact values 
fii oc and || /^ioc || 2 1 showing that suboptimal bounds still lead to meaningful bounds on the embedding 
dimension. 



5.2 



is a direct 



Remark 5.4. For fi < KN_^ 2 (as in Corollary 4.6), one has ||/i Zoc || 2 < K , so Theorem 
generalization of Corollary 4.6 As one has equality if and only if fi loc is constant, however, Theorem 5.2 
will be stronger in most cases. 



6 Local coherence estimates for frequencies and wavelets 

Due to the tensor product structure of both of these bases, the two-dimensional local coherence of the 
two-dimensional Fourier basis with respect to bivariate Haar wavelets will follow by first bounding the 
local coherence of the one-dimensional Fourier basis with respect to the set of univariate building block 
functions of the bivariate Haar basis. 

Lemma 6.1. Fix N — 2 P with p G N. For the space C N , the one- dimensional Fourier basis vectors ipk, 
k 7^ 0, and the one- dimensional Haar wavelet basis building blocks h e n kl e = 0, 1, satisfy the incoherence 
relation 

\(<p k ,h e nit )\ <min(^-^,37r2-5). 
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Proof. We estimate 



2P~ n gj r 2P~ n ~ 1 — 1 2 p ~ n £-\-2 p ~ n — 1 

2 p_„_l_ 1 



_ e 27ri2 (— l) e e 27ri2 ™ J 2*~ p - 

To estimate this expression, we note that 

|1 - e 2 * i2 ~ n ~ lk \ < min(2, vr2- n |fc|) (6.1) 

and distinguish two cases: 

If ^ |fc| < 2 P ~ 2 , we bound |1 - e 27ri2 ~ Pfc | > 2~ p \k\ and apply Q to obtain 

4 • 2t 

<min(— — ,27r2-i). 

For 2 p ~ 2 < |fc| < 2 P_1 , and hence 2~ p < we note that |1 — e 2?ri2 P,£ | > ^ an d bound, again using 



n „1 - e 



\ — g27ri2-Pfc 



2 



• /6- 22 „\ 

<mm( 3.2^). 



□ 

This lemma enables us to derive the following incoherence estimates for the bivariate case. 

Theorem 6.2. Let N — 2 P for N 9 p > 8. TTien £/ie focaZ coherence fi loc of the orthonormal two- 
dimensional Fourier basis {v?fei,fe 2 } wii/i respect to the orthonormal bivariate Haar wavelet basis {h e n e } in 
C NxN , as defined in (2.3) and (|2.2| 7 respectively, is bounded by 



M& 2 < K(h,k 2 ) := min (l, 

18ttV2 



< k'Qcx, k^) '■= min 1, 



12 i \U„\2\ 1 /' 2 



*1 3 + *2 



and one /ias ||ft||2 < ll^'lb < 52y/p = 52^/log 2 (-/V). 



Proof. First note that the bivariate Fourier coefficients decompose into the product of univariate Fourier 
coefficients: 

For fc, ^ 0, the factors can be bounded using Lemma [6. 1| which, for fci ^ 7^ yields the bound 

'6-2t n „\ . /6-2? „ „_»\ 18tt 



1(^,^,^)1 < min (^-^,3.2-t) min (^-^,3.2-t) < 



max(|fei|, |fc 2 |)' 



Next we consider the case where either k\ — or fc 2 = 0; w.l.o.g., assume ki = 0. We use that in one 
dimension, we have (<po, h x n ^ =0 as well as (ipo, /i° e ) = 2~5 . So we only need to consider the case that 
ei = and hence e 2 = 1. Thus we obtain 

., „ J'25 6 

KW),fa,^)| <2— 



|/c 2 | maxdfcilJfcal)' 
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In both cases, we obtain H l °^ k2 < max Q 1 ^ n^ry . The bound /^ fc ° c fc9 < 1 follows directly from the Cauchy- 
Schwartz inequality. We conclude fi l k oc k2 < n{kx,k 2 ) < n'(ki,k 2 ). 
For the Abound, we use an integral estimate, 



np— 1 

2 M8jT 2 



<#{(fc 1 ,fc 2 ):fc 1 +^< 648^}+ ^ , |2 + , 



|A-'ir + |fc2| 2 >6487r 2 



< 20600 + J J 187r\/2r _1 dr# 



r=187iV2-l 

< 17200 + 5021og 2 (AT) < 2700 log 2 (N) = 2700p, 

where we used that p > 8. Taking square root implies the result. □ 

As the infimum of a strictly decreasing function and a strictly increasing function is bounded uniformly 
by its value at the intersection point of the two functions, Lemma |6.1| also gives frequency-dependent 
bounds for the local coherence between frequencies and wavelets in the univariate setting. 

Corollary 6.3. Fix N = 2 P with p G N. For the space C N , the one- dimensional Fourier basis vectors 
ifk, k ^ 0, and the one- dimensional Haar wavelets satisfy the incoherence relation 

\(<Pk,K,i)\ < 3V2n/Vk. 



7 Recovery guarantees 
7.1 Proof of Theorem IQ1 



The proof of Theorem 3.2 concerning recovery from £i-minimization in the bivariate Haar transform 



domain follows by combining the local incoherence estimate of Theorem 6.2 with the local coherence 
based reconstruction guarantees of Theorem |5.2| Under the conditions of Theorem |5.2| the stated recovery 
results follow directly from Theorem |4.2| The weighted ^ 2 -norm in the noise model is a consequence of 
the preconditioning. 



7.2 Preliminary lemmas for the proof of Theorem 3.1 



The proof of Theorem |3. 1| proceeds along similar lines to that of Theorem 3.2 but we need a few more 
preliminary results relating the bivariate Haar transform to the gradient transform. Each of the following 
results are derived from a more general statement involving the continuous bivariate Haar system and 
the bounded variation seminorm. 

Proposition 7.1. Suppose f G , and suppose its bivariate Haar transform w — Hf G C N is 
arranged such that W(k) is the k-th largest-magnitude coefficient. Then there is a universal constant 
C > such that for all k > 1, 

\w(k) \ < C - 

In words, this proposition says that the bivariate Haar coefficient sequence of a function / is in weak 
i\ and its weak t\ semi-norm is bounded by the total variation semi-norm of /. Sec |25| for a derivation 



of Proposition 7.1 from Theorem 8.1 of [13] . 

We also have the following result about the bivariate Haar system. 

Lemma 7.2. Let N = TP . For any indices (ti,i 2 ) and (£i,£ 2 + 1), there are at most 6p bivariate Haar 
wavelets h e n( satisfying \h^ e (ti : t 2 + 1) - ^.^(^1,^2)! > 0. 

Proof. The lemma follows by showing that for fixed dyadic scale n in < n < p, there are at most 6 
Haar wavelets with edge length 2 P ~™ satisfying \h^ t {t\,t2 + 1) — /i^j f (ii,i 2 )| > 0. If the edge between 
(ti,t 2 ) and (ti } t 2 + 1) coincides with a dyadic edge at scale n, then the 3 wavelets supported on each of 
the bordering dyadic squares transition from being zero to nonzero along this edge. On the other hand, 
if (ti,t 2 ) coincides with a dyadic edge at dyadic scale n + 1 but does not coincide with a dyadic edge at 
scale n, then 2 of the 3 wavelets supported on the dyadic square having (ti,t 2 + 1), {ti,t 2 ) as a center 
edge can satisfy the stated bound. □ 
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Lemma 7.3. 

||V/i^||i<8 Vn,£,e. 

Proof. h e n e is supported on a dyadic square of side- length 2 p ~ n , and on its support, its absolute value is 
constant, \h e n g \ = 2 n ~ p . Thus at the four boundary edges of the square, there is a jump of 2 ra_p , at the 
(at most two) dyadic edges in the middle of the square where the sign changes there is a jump of 2 • 2 n ~ p . 
Hence HV/i^Ji < ||V/i^ 1} ||i < 8 • "2P~ n ■ 2 n ~ p = 8. □ 

We are now ready to prove Theorem |3.1| 
7.3 Proof of Theorem |3TH 

Recall that % : C N — > C N denotes the bivariate Haar transformation / >-> (^/i^^) n<E » let 
denote the j-th largest-magnitude Haar coefficient, and let denote the associated Haar wavelet. 

Let D e C N xN be the diagon al m atrix encoding the weights in the noise model, i.e., D 
diag(p), where, for k! as in Theorem 

|fc 2 | 2 ) 1/2 /187rY Note that Dg = pog. 



G.2 



p(k 1 ,k 2 ) = Wn'h/K'ikuki) = C v /log 2 (A r )max(l,(|fc 1 | 2 + 



By Theorem 5.2 combined with the bivariate incoherence estimates from Theorem 6.2 we know that 
with high probability A := -^^DTq/H.* has the restricted isometry property of order s and level i5 once 

m > s<5~ 2 log 3 (s)log 5 N. 

Thus, for the stated number of measurements m with an appropriate hidden constant, we can assume 
that A has the restricted isometry property of order 

s = 24<7 2 slog 3 (iV), 

where the exact value of the constant C will be determined below. In the remainder of the proof we show 
that this event implies the result. 



Let u = / — denote the residual error of (3.1 ). Then we have 

• Cone Constraint on Vu. Let S denote the support of the best s-sparse approximation to V/. 
Since /# = / — u is the minimizer of (TV) and / is also a feasible solution, 

ll(V/)s||i - IKVuJslli - ||(V/) s .||i + \\{Vu)sA\x 

< ||(V/) S - (Vu) s \U + ||(V/)s« - (Vu)s4i 

= I|V/ # || 1 

<\\Vf\\i 

= ll(V/) S ||l + ||(V/) S e|| 1 

Rearranging yields the cone constraint 

IKWMIi < ll(Vw)s||i + 2||V/ - (Vf)s\\i. (7.1) 



• Cone Constraint on w u = Hu. Proposition |7.1| allows us to pass from a cone constraint on 
the gradient to a cone constraint on the Haar transform. More specifically, we obtain 



7.2 



the set A index- 



Now consider the set S consisting of the s edges indexed by S. By Lemma 
ing those wavelets which change sign across edges in S has cardinality at most |A| = 6slog(A^). 
Decompose u as 

u = J2 w U) h U) =J2 w U) h U) + 12 w U) h U) = : UA + MAc 
j jeA jeA c 
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and note that by linearity of the gradient, 

Vll = VttA + V«A C - 

Now, by construction of the set A, we have that (Vua^s = and so (V«)s = (Vua)s- By Lemma 
|7.3| and the triangle inequality, 

||(Vu)s||i - ||(V« A )sr||i < ||Vua||i 

j£A 

<sEK-)l- 



Combined with Lemma 7.1 concerning the decay of the wavelet coefficients and the cone constraint 
(7.3), and letting 

s = 6slog(iV) = |A|, 
this gives rise to a cone constraint on the wavelet coefficients: 

N 2 N 2 

E H)\< E H)\ 

j=8+l j = s + l 

<Clog(iV 2 / s )||Vu||i 



= C\og(N'/ s )(\\(Vu) s \\ 1 + \\(Vu)s4 

< Clog(7V 2 / s )(||2(V U ) s ||i + 2||V/ - (V/) s ||i) 

< Clog(7V 2 /s)(l6E l w O)l + 2HW- (V/) s ||i 



J'SA 



< C\og(N 2 /s)^M + IIV/- (V/) s ||i 



3 = 1 



• Tube constraint, ||^l'Hu||2 < \2e. 

By assumption, A — -+^DJ 7 q'H* : C N — > C™ 1 has the RIP of order s > s. Also by assumption, 
\\DTnf — Dy\\2 — \\p° (J-nf — y)\U < \/me, so / is a feasible solution to (|3.1[). 



Since both / and are in the feasible region of (3.1 ), we have for u = / — 



m|MWu||| - IpJhW^Hl = ||23J=h«||i 

< \\DFnf - Dy\\\ + \\DTnf* - Dy\\\ 

< 2me 2 . 

• Using the derived cone and tube constrain ts on Hu along with the assumed RIP bound on A, the 
proof is complete by applying Proposition 4.3 using 7 = Clog(iV 2 /s) < 2Clog(A^), k = 6s log N, 
and £ = Clog(7V 2 /s)|| V/ — (V/)s||i- In fact, this is where we need that the RIP order is s, to 
accomodate for the factors 7 and k. □ 

8 Summary and outlook 

We established reconstruction guarantees for variable-density discrete Fourier measurements in both the 
wavelet sparsity and gradient sparsity setup. Our results build on local coherence estimates between 
Fourier and wavelet bases. The resulting sampling strategies are specific to two-dimensional discrete 
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images, that is, N x N blocks of pixels. A priori, our results do not directly generalize to higher 
dimensional or univariate signal models. In particular, optimal sampling strategies as well as stable 
image recovery guarantees remain open for higher-dimensional signals. 

Variable density sampling in compressive imaging has often been justified as taking into account the 
tree-like sparsity structure of natural images in wavelet bases (e.g., in [39]). We note that our theory does 
not directly exploit this additional structure, and depends only on the local incoherence between Fourier 
and wavelet bases. We expect, however, that this additional structure can be used to derive sampling 
strategies with stronger reconstruction guarantees, as indicated by the suboptimality of the sampling 
densities predicted by our results in numerical simulations (Figure [T]) . 

All the recovery guarantees in this paper are uniform, that is, we seek measurement ensembles which 
allow for approximate reconstruction of all images. For non-uniform recovery guarantees, we expect that 
the number of measurements required in our main results can be reduced by several logarithmic factors 
by following a probabilistic and "RIP-less" approach [5J. 

It should also be noted that this paper does not address the important issue of errors arising from 
discretization of the image and Fourier measurements. In particular, as observed for example in [I], 
the use of discrete rather than continuous Fourier representations can be a significant source of error 
in compressed sensing. The authors of [T] propose to resolve this issue using uneven sections, that is, 
the number of discretization points in frequency is chosen to be larger than the number of discretization 
points in time. Nevertheless, the results in [T] are again just formulated for incoherent samples. Recently, 
it has been proposed to overcome this issue by sampling all of the low frequencies in addition to uniformly 
sampling the higher frequencies [2], but to date, no provable reconstruction guarantees have been provided. 
We expect that our approach can be applied to this setup - due to the variable density, it may even be 
possible to sample from the infinite set rather than restricting to a finite subset based on an intricate 
criterion. Such a generalization is out of reach for the optimization-based approaches such as in [30] . 
which will always be specific to the given problem dimension. In this sense, we expect that the additional 
understanding provided by this paper can eventually lead to optimized sampling schemes. All these 
questions, however, are left for future work. 



A Proof of Proposition 4.3 



Write use = u' 1 ) + + • • • + where r = Lip^J- Here vS l > is the 4fc7 2 -sparse image consisting of 
the 4&7 2 largest-magnitude components of us<=, consists of the 4fc7 2 largest-magnitude components 
among the remaining entries of us<=, and so on. Since the magnitude of each component of uC? — *) is at 
least as large as the average magnitude of the components of , we obtain 

11^0-1)11, 

Combining this with the cone constraint gives 



u) h < — ^IMIi < ^IMIi + ^7f^ < hush + iAtt^- 



Yuu , . .... . 



Together with the tube constraint and the RIP, we obtain 



> \\Au\\ 2 



> \\A(u s + u^)h-J2\\A(u^)h 



3=2 



> VT^5\\u s + u {1) \\ 2 - Vl + 5j2\\u^\\ 2 

3=2 

> VT^~s\\u s + u^\\ 2 - VTT~s(l\\u s \\ 2 + ^- / =^ 
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Then, since 6 < 1/3, 



||«s + u (1) || a < 5e + 



3£ 

"fVk 



Finally, because HE^^IIa < E5=a ll« W) lla < §K + ^ {1) h + we have 

||u|| 2 <fe+^L 
7V k 



confirming (4.3 ) 



□ 
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1G 



(a) 



(b) 




(g) 




Fi gur6 1; Various reconstructions of an MRJ image x £ ]]^256x256 w ^ n total variation minimization as in Theorem |3.l| 
with e = and using m = 6400 noiseless partial DFT measurements y = Fqx with frequencies Q = (fci,Ai2) sampled 
from various distributions. Beside each reconstruction is a plot of K-space {(fci,A;2) : —N/2 + 1 < fei,fe2 < N/2} and 
the frequencies used in its reconstruction, (a) Original image (and all of K-space). (b) Reconstruction using only lowest 
frequencies: Q = {(fci, fea) : fc^+fc| < 80}. (c) Prob((fci, k 2 ) S f2) ~ 1 (Uniform subsampling) (d) £1 comprised of frequencies 
in equispaced radial lines, (e) Prob((fci , fc 2 ) £ ft) oc (k\ + k\ + l) _1 / 2 (f) Prob((fci , fc 2 ) £ !1) oc (max(|fci|, |fc 2 |) + 
(g) Prob((fci, k 2 ) e n) oc (k'f + k\ + l) -1 . (h) Prob((fci, k 2 ) 6 fl) a {k'{ + k\ + 1)~ 3 / 2 . Theorem 3.1 guarantees stable 
and robust recovery for the inverse square-distance distribution in (g); a slightly stronger guarantee can be obtained for 
the inverse-max sampling distribution given in (f) from the stronger local coherence bound in Theorem |6.2| The relative 
reconstruction error ||/ — /*y||2/||/||2 corresponding to each reconstruction is (b) .2932, (c) .8229, (d) .4074, (e) .3192, (f) 
.2603, (g) .2537, and (h) .2463. 
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Figure 2: The same reconstructions as in Figure [l] magnified to the chin region: (a) Original image (and all of K- 
space). (b) Reconstruction using only lowest frequencies: Q = {(ki,k2) : k\ + fc| < 80}. (c) C comprised of frequencies 
in equispaced radial lines, (d) Prob((fci , fc 2 ) S £1) oc (kj + k\ + l) -1 / 2 (e) Prob((fci , fc 2 ) £ fi) tx (fc 2 + fcf + (f) 
Prob((fei,fc 2 ) eSJ)« (fcf + k\ + l)" 3 / 2 . Theorem [3T 



concerns stable and robust recovery for method (c). 
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